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Abstract 

Recent achievements in experiments with cold fermionic atoms indicate the poten- 
tial for developing novel superconducting devices which may be operated in a wide 
range of regimes, at a level of precision previously not available. Unlike traditional, 
solid-state superconducting devices, the cold-atom systems allow the fast switching 
on of the BCS phase, and the observation of non-equilibrium, coherent oscillations 
of the order parameter. The integrable and non-linear nature of the equations of mo- 
tions makes this operating regime particularly rich in potential applications, such as 
quantum modulation and encoding, or nonlinear mixing of quantum coherent oscil- 
lations, to name only two. From a mathematical point of view, such systems can be 
described using the Knizhnik-Zamolodchikov-Bernard equation, or more generally, by 
the mattix Kadomtsev-Petviashvilii integrable hierarchy. This identification is partic- 
ularly useful, since it allows a direct application of the known non-linear phenomena 
described by particular solutions of these equations. Other important features of this 
formulation, such as the relation to the spin Calogero-Sutherland model, also have 
relevant physical interpretations. In this work, a complete description of these rela- 
tionships is presented, along with their potential practical consequences. 
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1 Introduction 

The operating parameters specific to experiments with cold fermionic gases exhibit- 
ing the paired BCS state [31-38] allow the investigation of previously inaccessible 
regimes, as suggested in [17]. It was shown that for such systems, the time scales of 
the order parameter ta ~ | A | ~ 1 , and the quasiparticle energy relaxation time r c are 
both much larger than typical time for switching on the pairing interaction r , essen- 
tially given by the variation of external parameters, such as detuning from the Feshbach 
resonance. It was argued that in this regime, for times t <C r e , the dynamics of the 
system is given by non-linear, non-dissipative equations describing the coherent BCS 
fluctuations for the system out of equilibrium. In this limit, the system is integrable, 
and features non-perturbative behavior, such as soliton-type solutions. 

In the mean-field limit, such non-trivial solutions describing the collective mode 
of the Anderson spins [39] were derived in [17], for a two-level effective system. This 
work was generalized [15] in algebro-geometrical terms. 

In [40—42], the long-time behavior of the solution has been considered, under vari- 
ous conditions. An issue that remains incompletely addressed so far is the relaxation of 
the nonlinear oscillatory solution induced by perturbations of the spectral curve, phys- 
ically justified by coupling to the environment. Several possible kinds of perturbations 
may be considered, which may lead to different types of relaxation. 

In this article, we describe the mathematical structure underlying the integrable 
quantum pairing problem, and its relations with the theory of isomonodromic defor- 
mations of differential equations and other integrable models of strongly interacting 
particles. We also explain the classical limit of the quantum hamiltonian, and give 
the explicit construction of the solution in this case. Possible applications due to the 
integrability of the model are listed in the conclusion of the article. 

The paper is organized as follows: in section|2j the Richardson-Gaudin model and 
its relation to other integrable models are reviewed. Section [3] describes the classical 
limit and explicit multi-phase solution. In section H] a few possible applications are 
briefly discussed. 
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2 The Richardson-Gaudin Model 

As argued for the first time in [17], optically trapped cold Fermi gases in the BCS 
state are correctly described by the Richardson-Gaudin quantum pairing model. Using 
this model is justified by the discreteness of the energy spectrum, and by the spatially 
uniform profile characterizing such systems. This integrable system, introduced by 
Richardson and Sherman [1-6] in the context of nuclear physics, has received revived 
interest in recent years, after being applied to metallic superconducting grains [16]. 
The model is intimately related [18-20] to a class of integrable systems generally 
referred to as Gaudin magnets [7]. These systems have been studied both at quantum 
and classical level [21-30], in the elliptic case as well as trigonometric and rational 
degenerations, using various methods from integrable vertex models to singular limits 
of Chern-Simons theory. 

In this section, we recall the main features of this model, as well as their physical 
interpretations. 



2.1 The quantum pairing hamiltonian 

Following [21], we briefly review the Richardson pairing model. It describes a system 
of n fermions characterized by a set of independent one-particle states of energies e ; , 
where the label I takes values from a set A. The labels may refer, for instance, to 
orbital angular momentum eigenstates. Each state I has a total degeneracy di, and the 
states within the subspace corresponding to I are further labeled by an internal quantum 
number s. For instance, if the quantum number I labels orbital momentum eigenstates, 
then di = 21 + 1 and ,s = —I, . . . ,1. However, the internal degrees of freedom can be 
defined independently of I. We will assume that di is even, so for every state (Is), there 
is another one related by time reversal symmetry (Is). For simplicity, we specialize to 
the case di — 2, s =| , j. Let c] represent the fermionic creation operator for the state 
(Is). Using the Anderson pseudo-spin operators [39] (quadratic pairing operators), 
satisfying the su(2) algebra 
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the Richardson pairing hamiltonian is given by 



Hp = 2 ^f - a J2 = E 2eit i - a t+ ■ *" 



(1) 



(2) 



i,v 



leA 



where t = tj is the total spin operator. It maps to the reduced BCS model 
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by replacing the translational degrees of freedom by rotational ones, where I 6 A = 
{1, . . .n} ennumerates the one-particle orbital degrees of freedom, while s = j, J, 
indicates the two internal spin states per orbital (di = 2). The pairing hamiltonian can 
be decomposed into the linear combination 
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At a fixed value of the component t 3 of the total angular momentum, the last term be- 
comes a constant and is dropped from the hamiltonian. The operators Ri (generalized 
Gaudin magnets [7]) are given by 

These operators solve the Richardson pairing hamiltonian because [18] they are inde- 
pendent, commute with each other, and span all the degrees of freedom of the system. 
Richardson showed [1,2] that the exact N— pair wavefunction of his hamiltonian is 

t t 1 

given by application of operators b\ = J2 t 2ei -e k to vacuum (zero pairs state). The 
unnormalized N— pair wavefunction reads '^(e;) = Yik=i fyt|0)- The eigenvalues e% 
satisfy the self-consistent algebraic equations 



which can be given a 2D electrostatic interpretation [21] with energy 
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Equations © appear as equilibrium conditions for a set of charges of strength q = 2 
placed at points efe, in the presence of fixed charges of strength q = —1 at points 2e/, 
and uniform electric field of strength i, pointing along the real axis. This interpreta- 
tion proves to be very useful for the conformal field theory (CFT) description of the 
Richardson problem. The electrostatic energy (0 is minimized for values {e^} corre- 
sponding to pair energies. In (0, n, N represent the number of single-particle levels 
and the number of pair energies, respectively. For large interaction constant g, the 
equilibrium positions {et} form a set of complex conjugated pairs defining a curve 7 
in the complex plane of energies. We note that the eigenvalues r» of Gaudin hamilto- 
nia are proportional in this language to the values of the electric field at positions 2ej, 

For a set of single-particle energies {ei}, the BCS ground state is obtained by 
minimizing the electrostatic energy with respect to positions of the free charges 
at {efc}. Once found, they also determine exactly the values of the electric field at 
positions {2e, } on the real axis, which are proportional to the ground-state eigenvalues 
{ r i }• F° r any other values of {r^}, the electrostatic energy (0 is not minimized. This 
indicates that for arbitrary values r* ^ rf s , the system is not in equilibrium. 



2.2 Richardson- Gaudin model as limit of KZB equations 

In this section we recall the relation between the Richardson-Gaudin model and sin- 
gular SU(2) Chern-Simons and WZW models [25,26]. These relations stem from the 
study of quantum states of the Chern-Simons theory with gauge group SU (2) on the 
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manifold T 2 x R and in the presence of Wilson lines {zt} X M, i = 1, . . . , n [25], 
corresponding to the one-particle spectrum of the Richardson hamiltonian, Zi = 2e;. 
The torus T 2 has modular parameter r = t\ + iT2, rj.,2 G JR. Thus, T 2 = C/Z + rZ. 
The quantum states of this theory are known to satisfy the Knizhnik-Zamolodchikov- 
Bernard (KZB) equations 

W C s = 0, V = (V T ,V Zi ), (9) 
where V is the flat KZB connection, 

V T = nd T + H (t, Zi ), V Zi = nd Zt + Hi(r, z,), (10) 

and Ho , Hi are elliptic versions of the Gaudin hamiltonia. The parameter n is related 
to the level of representation of SU(2). By a limiting procedure, the system degen- 
erates into the Richardson-Gaudin problem: by taking the limit r — > ioo, the torus 
degenerates into a cylinder, and the elliptic hamiltonia Hi,Ho take a trigonometric 
limit. Upon rescaling the cylinder such that the set {zi} is collapsed into the origin, 
the cylinder becomes the complex plane with punctures at {z^}, and the rescaled op- 
erators become the rational Gaudin magnets ([5j. 

2.3 Singular representations of Wess-Zumino-Witten theory 

Given the well-known relation between 2 + 1 Chern-Simons theory with SU (2) gauge 
group and WZW theory, it is not surprising that a non-Abelian CFT representation for 
Richardson-Gaudin models is also available. In fact, the relationship between these 
theories is already transparent by identifying the electrostatic picture of R-G with a 
Coulomb gas representation of WZW. 

In the CFT approach [16], a primary field of spin 1/2 of the SU(2) representation 
of WZW is associated to every energy level of the R-G problem. The background 
charge a and the level of Kac-Moody algebra SU (2)*. are related by k = tt-^t — 2. It 
turns out that the R-G model is retrieved in the singular limit a a — > oo, corresponding 
to k = —2. In this limit, the representations are given by perturbed WZW blocks 
(PWZW), 

* P (zi)= <f de l( f de N e- a o u ^ ek ^ R { Zi ,e k ). (11) 

J Ci ^ Cjv 

In ( fTTl i. efc are the eigenvalues of the Richarsdon hamiltonian with N fermion pairs, 
ao is the background charge in the Coulomb gas representation of the WZW model 
associated with the 2+1 Chern-Simons theory described above, and U(zi, e^) is the 
holomorphic piece of the 2D electrostatic potential of the Coulomb gas representation 
(0. The Richardson wave function ^n(zi,ek) is obtained by applying Bose fields 
corresponding to positions Zi,ek, 

n N 

^R(z i ,e k ) = (l[^(z i )l[(3(e k )), (12) 
i=l fc=l 

and with correlator (P(z)-f(z')) = jz^j- 

The singular limit described above in the language of the Chern-Simons theory 
corresponds to the limit ao — > oo in the CFT approach. The perturbed WZW blocks 
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are dominated by the saddle-point configurations, given by the eigenvalue equations 
d ek U = 0, and become proportional to the Richardson wavefunction. In the singular 
limit, the KZB equations takes the form 



which shows that the Richardson wavefunction is an eigenstate of the Gaudin hamil- 
tonia, with eigenvalues given by the solutions to the electrostatics problem (0. 

Equation ([T3l is the analog of the singular limit of the KZB equation (0, in the 
limit k — ► 0, corresponding to a — > oo. Likewise, there is a dirrect correspondence 
between the Chern-Simons states and the perturbed WZW blocks. We will see in the 
following sections that these equations survive in the semiclassical limit and that they 
allow the identification of the tau function for the corresponding classical integrable 
system. 

2.4 Relation with the KP hierarchy and other integrable models 

The zero curvature conditions for the KZB connection can be interpreted as a trun- 
cation of the general Kadomtsev-Petviashvilii (KP) integrable hierarchy, at a lelvel 
set by the number of eigenvalues in the single-particle spectrum. This is also illus- 
trated by the fact that, as shown in the preceding paragraphs, in the genus one case, 
the Richardson-Gaudin problem is equivalent with the elliptic Calogero model. This 
model is known to be embedded into the KP hierarchy, and admits simple-pole solu- 
tions built from the solution of KP. In fact, all the relationships described so far are 
consequences of an even more general formulation of the Richardson-Gaudin problem 
as a hierarchy of isomonodromic deformations with respect to changing the moduli of 
its spectral curve. 

These relationships allow to classify the solutions to the quantum Richardson- 
Gaudin model using degenerate limits of the solutions for the elliptic Calogero prob- 
lem. Such solutions are known (for instance, the simple-pole solutions mentioned 
above) and would therefore allow (in principle) constructing exact solutions for given 
initial conditions. Even more importantly, there is a large body of knowledge con- 
cerning the classes of solutions supported by the KP hierarchy (periodic, algebraic- 
geometric or rational), which translates into distinct "phases" (albeit we are dealing 
with non-equilibrium, coherent oscillations) supported by the quantum pairing model, 
depending on the initial conditions. 

3 The mean-field limit of Richardson-Gaudin models 

Throughout this section, we consider the classical limit of the quantum pairing model. 
The main reason for this approximation is the fact that explicit, generic solutions can 
be obtained in this case. 

3.1 General description of the classical model 

In the mean-field limit, the spin operators t; are replaced by their quantum mechan- 
ical averages. Written in terms of the classical vectors S; = 2(t;}, the semiclassical 
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Figure 1: Schematic representation of the Liouville torus for the integrable system ([TBI . 



approximation for the pairing hamiltonian becomes 
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where J = J2ie\ ^ l m ^ tne ^CS gap function is given by A = gj /2. Replacing 
commutators by canonical Poisson brackets, 



{Sf, Sf] = 2e a ^S]5 ih 



(15) 



variables Sf become smooth functions of time. In this limit, the problem can analyzed 
with tools of classical integrable systems, and the solution is known to be exact as 

n — > oo. 

The Poisson brakets ( fT3] l and hamiltonian (TBI) lead to the equations of motion 



S, = 2(-A + e t z) x Si, 



(16) 



where 2 A = (gj x ,gj y , 0) and J is the total spin. The semiclassical limit of Gaudin 
hamiltonia are independent constants of motion, 



1 

T'l = — 

2 
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Si • Sj 



h = o. 



(17) 



Equations ( fTBT ) describe a set of strongly interacting spins and have generic non-linear 
oscillatory solutions. The exact solution may be obtained through the Abel-Jacobi 
inverse map, explained in the next section. 



3.2 Isospectral case and the Abel-Jacobi inversion problem 

In [9, 15], the system (fT~6b with fixed spectral curve was solved through inverse Abel- 
Jacobi mapping, by using Sklyanin separation of variables techniques [8,52]. Interest- 
ing connections to generalized Neumann systems and Hitchin systems were discovered 
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in [24]. The solution starts from the Lax operator 



C(X) = -a 3 
9 



E 



Si • ' 
A^ 



a(A) 6(A) 
c(A) -a(A) 



(18) 



where er Q , a = 1, 2, 3 are the Pauli matrices, and A is an additional complex variable, 
the spectral parameter. Let Uk, k = 1, . . . , n — 1 be the roots of the coefficient c(A). 
Poisson brackets for variables Sf read 



{Sf, Sf } = 2e^ 7 ^5 



(19) 



The Lax operator ( fT8l defines a Riemann surface (the spectral curve) T(y, A) of genus 
g = n — 1, through 



P(A) 



n 2 



?/ = Q(A) =det£(A) 



whereP(A)=nr= 1 (A-e i )- 

The equations of motion for the hamiltonian ( fT4b become 



(20) 
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u- = r 



9 J 3 + 2J2 



fe=i 



(21) 



In (ED, u = -2 EL - ! 1 u h b( Ui ) = 0. 

From the equations of motion, it is clear that knowledge of the initial amplitude 
of J~ and of the roots {ui} is enough to specify the n unit vectors {S^}, for a given 
set of constants of motion {Ri} given by the classical limit of Gaudin hamiltonia. 
The Dubrovin equations (f2Tb are solved by the inverse of the Abel-Jacobi map, as we 
explain in the following. We begin by noting that the polynomial Q(\) has degree 
2n, and is positively defined on the real A axis. Therefore, the curve T(y, A) has n 
cuts between the pairs of complex roots [i^-i, E2i], i — 1,2, ... ,n, perpendicular 
to the real A axis. The points u; belong to n — 1 of these cuts, ui G [Ezi-i, E2i], i = 
1, . . . , n — 1. These g — n — 1 cuts allow to define a canonical homology basis of 
r, consisting of cycles {cti, fa}, i = 1, . . . ,g. With respect to these cycles, a basis of 
normalized holomorphic differentials {oji] can be defined, through 



Mi = A 9_i — , My = 
V 



Hi, u> = M jLt. 



(22) 



The period matrix By — j^ .oJi is symmetric and has positively defined imaginary 
part. The Riemann 8 function is defined with the help of the period matrix as 



(z\B) = ^2 e 27Ti{ntz+ ^ ntBn) . 



(23) 
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The g vectors Bk consisting of columns of B and the basic vectors define a lat- 
tice in C 9 . The Jacobian variety of the curve T, is then the g— dimensional torus 
defined as the quotien J(T) = C 9 /(Z 9 + BZ, 9 ). The Abel-Jacobi map associates to 
any point P on T, a point (g— dimensional complex vector) on the Jacobian variety, 
through A(P) = J P lj. Considering now a g— dimensional complex vector of points 
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{Pfe}, k = 1, . . . , g on r, defined up to a permutation, we can associate to it the point 
on the Jacobian 

a 

z = a(P)=J2 A (Pk)+K, (24) 

fe=i 

where K is the Riemann characteristic vector for V. 

The map d24b suggests that we now have a way to describe the dynamics on T 
by following the image point on the Jacobian. Given a point on the g— dimensional 
Jacobian z = . . . , Cn-iX we can find an unique set of points {Xk}, k = 1, . . . , g 
on r, such that z = a(\), and 9(a(P) — z\B) = 0. The system evolves in time 
according to the point z(t) 

Ck = ick, I < i < g - 1, Cn-i = i( c n-i + t), (25) 

where {ck} is a set of initial conditions, such that Zq = z(t = 0) = a(c), and 
c is the set of initial conditions for positions of A on T. Together with the initial 
condition which determines the initial amplitude of J~, this set will determine entirely 
the evolution of the functions u,i{t) 1 J~ (t). 

It is perhaps useful to give a graphical description (FigureQ]) of the general solution 
described above. For a given set of initial conditions for the spins {Si}, i.e. also of 
the constants of motion {r.;}, a polynomial Q(u) of degree 2n and with n pairs of 
complex conjugated roots E 2k+2 = E 2 k+i, k = 0, ... ,n — 1, is constructed. A 
schematic representation of these roots is given in Figure Q] Between each pair of 
roots, we place a simple cut Ck = [E^fc+i, E 2 k+2] on the complex plane of energies. 
The surface thus obtained is a representation of a torus of smooth genus g = n — 1. 

Variables Uk are introduced for re — 1 of these cuts, with respect to which the 
equations of motion separate. The variables Uk evolve in time in a complicated fashion, 
solving a system of nonlinear coupled differential equations ( f2TT >. Up to a constant, the 
time dependence of the gap parameter amplitude is given by 

log|A(£)| =lm / u(t)dt, u(t) =^u fc (i). (26) 
^ k 

The widths of the cuts Ck and the periods of the non-linear oscillators Uk are deter- 
mined by the values of constants of motion {r;}. For the particular choice n — rf s , 
all the cuts Ck, k = 1, . . . , n — 1 vanish, and the width of the remaining cut C n equals 
the equilibrium value of the gap function: \E 2 n~i — E 2 ^\ = 2|A| GS . In that case, 
the oscillators Uk = ^2fe-i = E 2 k are at rest, and the only time dependence left in 
the system is the uniform precession of the parallel planar spins £r, with frequency 

ijj = 2 J2k=i € k — 2 J2pZi E 2p -i — Y17=i In the case of particle-hole symmetty, 
w vanishes as well. 

4 Conclusion 

In the present work, the regime of coherent, non-equilibrium quantum oscillations 
characterizing the onset of the BCS state in cold Fermi gases was investigated from a 
mathematical point of view. The integrability of the model and main features of the 
space of solutions were described, both at the quantum and classical levels. 
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The most direct application of these coherent quantum oscillations would be in 
encoding/decoding of information. Although simple in principle, this would in fact re- 
quire very accurate control over the initial conditions at the onset of the pairing regime, 
which is unfeasible. However, given the stability of the coherent oscillatory regime, it 
is clearly easier to first implement the quantum modulation of an arbitrary oscillatory 
solution, by slowly perturbing the single-particle spectrum. This kind of perturbation 
is in fact within the same class as the original isomonodromic deformations, and hence 
does not destroy the integrability of the model. 

Another possible application makes use of the non-linear nature of the solution 
of the Richardson-Gaudin model. Considering only the simplest classical solutions 
given by the Abel-Jacobi construction, the elliptical solutions, it is possible to devise 
a non-linear 'mixer' by using the fact that the dynamics is linear on the associated 
torus (Jacobian). Then, linear operations with respect to the dynamics on the Jacobian 
would translate into strongly non-linear (but still invertible) effects observed at the 
level of quantum oscillations. 
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